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Abstract 

Spectral integration is a method for solving linear boundary value problems which uses the Chebyshev 
series representation of functions to avoid the numerical discretization of derivatives. It is occasionally 

T-H attributed to Zebib (J. of Computational Physics vol. 53 (1984), P- 443-4-55) and more often to 

Greengard {SIAM J. on Numerical Analysis vol. 28 (1991), p. 1071-1080). Its advantage is beheved 
to be its relative immunity to errors that arise when nearby grid points are used to approximate 
derivatives. In this paper, we reformulate the method of spectral integration by changing it in four 

h^- different ways. The changes consist of a more convenient integral formulation, a different way to treat 

and interpret boundary conditions, treatment of higher order problems in factored form, and the use of 

i-S^ piecewise Chebyshev grid points. Our formulation of spectral integration is more flexible and powerful 

d as show by its ability to solve a problem that would otherwise take 8192 grid points using only 96 

grid points. Some of the intermediate quantities have large errors which cancel in the final result. The 

I— —I matrices that arise will have large condition numbers if the solution develops thin boundary layers. 

However, spectral integration is accurate in spite of the ill-conditioning of those matrices. 
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1^ In this introduction, we begin by describing the method of spectral integration as it is ordinarily 

understood followed by the four modifications and enhancements we make to the method. For 

CN expository convenience, we consider the second order boundary value problem 

T— I 

> D'^u + bDu + cu = f, n(±l) = (1.1) 

^ where D = f{y) is known over the interval — 1 < y < 1, and u{y), — 1 < y < 1, is the 

^ solution to be determined. The coefficients b and c are constant scalars. The solution to this 

problem can develop boundary layers. For example, when 6 = and c is a large negative 
number, the width of the boundary layer at y = ±1 is typically and approximately (— c)~^/^. 

The method of spectral integration as described by Zebib jl61 1984] and Greengard [H 
1991] begins by expanding D^u in a truncated Chebyshev series of the form 



— o — + 2^ + — o — • ^^-^1 



The Chebyshev polynomial Tj{y) is defined for — 1 < y < 1 as cos jO if cos = y and < ^ < vr. 
The M Chebyshev grid points are cos(j7r/M), j = 0, . . . , M. The discrete cosine transform 
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is an efficient method for passing back and forth between function values at Chebyshev grid 
points and the Chebyshev coefficients. For convenience, we use the notation Tj{f) for the 
j-th coefficient in the truncated Chebyshev series of /. Using this notation, we may write 
aj = Tj {D'^u). In implementation, it is often convenient to set au = 0. 
The indefinite integral of Tn{y) is 



!n+l 



2(n+l) 2(n-l) 



if n > 1, 



T2/4 if n = 1, and Ti if n = 0. Likewise, the indefinite integral iterated twice is 



tn+2 



Tn 



+ 



4(n + l)(n + 2) 2(n2 - 1) 4(n-l)(n-2) 



if n > 2, 



(1.3) 



(1.4) 



r2/4if n = 0, r3/24 if n = 1, and Ti/AS-T2/Q if n = 2. The Chebyshev series oi D'^u+h Du+cu 
is obtained using the series (1.2) oi D^u and the indefinite integrals (1.3) and (1.4). The series 
of D'^u + b Du + c It is in terms of the coefficients aj which appear in (1.2) and will have a term 
of the form A + By corresponding to constants of integration. If the truncated Chebyshev 
series of / is set equal to the Chebyshev series of D^u + b Du + c n, we get M + 1 equations 
for the M + 3 unknowns oq, . . . , cxm, A, and B. This set of M + 1 equations is banded. The 
boundary conditions give two more equations [3]: 




+ 



- 1) 2{p - 1) 4j(j + 1) 



+ 



4j(j - 1) 2{p - 1) + 1) 



1)^ =0 



(1.5) 



where a'- = aj for j = 0, . . . , M — 1, 



aA//2, and a'j^. 



0. After these two 



equations are inserted, we have M + 3 equations for M + 3 unknowns but the system is no 
longer banded. This completes our description of the method of Zebib [16] and Greengard 

The first change we make to the method is as follows. The method of Zebib and Greengard 
begins with the Chebyshev series of D^u. Instead, we begin by integrating the differential 



equation (1.1) to get 



u+b u+c u+A+By 



f u{±l)=0. 



[1.6) 



Here the integrals are understood to be indefinite integrals. If we expand n as a Chebyshev 
series of the form (1.2), instead of expanding D^u as in the method of Zebib and Greengard, 



we may again use standard formulas to calculate the coefficients in the Chebyshev series of the 



left and right hand sides of (1.6) and then equate them. The boundary conditions now take 
on the form 



+ 



M-l 



+ 



aM 



2 



M-l 



0. 



^1.7) 



As before, the resulting matrix system is banded except for the boundary rows, which are 
dense. The matrix system that results was presented on page 119 of Gottlieb and Orszag [HI 
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1977] . Gottlieb and Orszag noted that the matrix system was "essentially diagonally dominant" 
and therefore well conditioned, neither of which is in fact true if the solution develops thin 
boundary layers, but did not explain how they derived it. However, the method has been used 
extensively in fluid mechanics for decades [6l [7]. Muite [9] has recognized that the matrix 
system of Gottlieb and Orszag is obtained if differentiation matrices are used to discretize 



(1.1) and the resulting discretization is multiplied by an integration preconditioner of the type 



discussed by Coutsias et al. [T]. It is advantageous to begin by reformulating the boundary 



value problem in the integral form (1.6) for developing other aspects of spectral integration. 

The second change we introduce improves the way boundary conditions are handled and 
brings to light a numerical property of spectral integration which a good implementation must 
have. In integral form, the boundary value problem is 

u + bfu + cffu + A + By= Iff. (1.8) 



If the boundary conditions m(±1) = are used directly the matrix system loses its banded 



structure. Instead we find a particular solution of (1.8) with the global conditions 



To{u) = Tiiu) = 0. 

These conditions are global because the first two Chebyshev coefficients To and 7i are defined 
by integrating over the entire domain as in J^^ u{y){l — y'^)~^^'^ dy and J^^ yu{y){l — y^)"^''^, 
respectively. The advantage is that these integral conditions lead to a banded system. Two 
homogeneous solutions are found by solving 



u + b J u + c J J u + A + By = 

subject to the global conditions 

Toiu) = 1, Ti{u) = and %{u) = I, Ti{u) = 0, 

respectively. Each homogeneous solution is found by solving a banded system. The boundary 
conditions u(ibl) = are enforced by adding an appropriate linear combination of the two 
homogeneous solution to the particular solution. 

It may seem that two homogeneous solutions can be found much more easily using the roots 
of the characteristic equation of + 6A + c = 0. While finding homogenous solutions using 
the characteristic roots is easy, doing so will ruin the accuracy of the method. Typically the 



particular solution of (1.8) satisfying the global or integral conditions Tq{u) = 7i(n) = will 
have a Chebyshev series that converges much more slowly than the particular solution which 
satisfies boundary conditions such as u(±l) = 0. If the homogeneous solutions are found using 
the characteristic roots, one has to use enough grid points to resolve the particular solution 
satisfying the global conditions Tq{u) = Ti{u) = 0. On the other hand, if the homogeneous 
solutions are found as indicated, the grid need only resolve the solution that satisfies the 
given boundary conditions. We show in Section 4, that errors in the particular solution and 
homogeneous solutions found using global or integral conditions cancel each other when those 
solutions are combined to find a solution u that satisfies ^(±1) = (or some other boundary 
conditions of that type) . 
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Fig. 1.1: Solution of {D^ - aD)u = with a = 10^ and u{l) = 2 and u{-l) = 1. The figure 
shows the very thin boundary layer near y = 1. The computed solution (solid line) is 
in excellent agreement with the exact solution (filled markers). 

The third change we make is to tackle higher order problems by factorizing the operator. If a 
linear boundary value problem can be written in the form {D — ai){D—a2){D'^-\-biD-\-ci)u = /, 
for example, the method we describe in Section 3 tackles the boundary value problem by 
considering the operators (D — ai), (D — 02), and {D"^ + biD + ci) separately. This modification 
eases the burden of implementation greatly. In principle, the generalization of the method of 
spectral integration to higher order problems is not difficult as we show in Section 2. However, 
the implementation becomes laborious. The factored version of spectral integration is described 
in Section 3. The idea of factorizing operators occurs in the work of Kim, Moin, and Moser 
j6], a landmark in the modern development of fluid mechanics, in a special form. It is used in 
a more implicit way in the earlier and equally important work of Kleiser and Schumann |7j . 

The fourth change, which is described in Section 7, is perhaps the most significant of all. 
In Section 7, the factored form of spectral integration is generalized to piecewise Chebyshev 
grids. Figure |1.1| shows a solution with a thin boundary layer. It takes 8192 grid points 
to compute the solution with 10 or more digits of accuracy using any of the usual forms of 
spectral integration. Spectral integration with piecewise Chebyshev grid points reaches 10 or 
more digits of accuracy using only 96 points. 

The example used in Figure |1.1| is taken from an important paper by Greengard and 
Rokhlin |5j (also see Starr and Rokhlin [12] )• In that paper, Greengard and Rokhlin show 
how to make numerical methods based on Green's functions work for problems with non- 
constant coefficients. For a boundary value problem of the form u" + p{y)u' + <l{y)u = f{y) 
with boundary conditions aX y = a and y = c, Greengard and Rokhlin split the interval 
[a, c] into subintervals each of which is discretized using Chebyshev points. The solution u 
is represented using a background Green's function. The method of Greengard, Rokhlin, and 
Starr appears to be the best general purpose method for solving linear boundary value problems 
in one spatial dimension. It can handle non-constant coefficients, boundary layers, and rapid 
oscillations. Spectral integration with piecewise Chebyshev grid described in Section 7 requires 
the coefficients to be constant, but is much cheaper when coefficients are in fact constant, as 
they are in many applications. In particular, spectral integration with piecewise Chebyshev 
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grid is likely to be of much importance for high Reynolds number channel and Couette flow 



When the boundary value problem admits thin boundary layers, the matrices that arise in 
spectral integration are very ill-conditioned. However, when the width of the boundary layer 
is 10~® and the condition number is as high as 10^^, only six digits of accuracy are lost. We 
investigate this partial immunity to poor condition numbers in Section 5. 

The numerical examples in Section 6 illustrate three properties of spectral integration: 
ability to resolve boundary layers, impressive accuracy even when the matrices that arise are 
ill-conditioned, and robustness when the solution is highly over-resolved. 

In addition to the generalization of spectral integration to piecewise Chebyshev grids. Sec- 
tion 7 derives a method for solving boundary value problems using piecewise Chebyshev grids 
and differentiation matrices. Spectral methods are much easier to implement when differenti- 
ation matrices are used explicitly. The disadvantages of using differentiation matrices are high 
cost, the number of arithmetic operations being cubic in the number of grid points as opposed 
to the linear cost of spectral integration, and large rounding errors when the number of grid 
points is increased. However, the use of piecewise Chebyshev grids appears to eliminate the 
need for a Chebyshev grid with a lot of points in many if not most problems. In Section 7, 
we show that spectral differentiation matrices are effective for solving problems with boundary 
layers and internal layers. The cost remains higher than that of spectral integration but the 
rounding errors are no longer a problem. 

2 Spectral integration of boundary value problems 

Section 2.1 describes spectral integration of first and second order linear boundary value prob- 
lems. The factored form of spectral integration described in Section 3 reduces higher order 
boundary value problems to such first and second order problems. Section 2.2 describes r-th 
order spectral integration, which tackles r-th order problems directly, and Section 2.3 is a 
reformulation of the method of Zebib and Greengard. Both r-th order spectral integration 
and our formulation of Zebib-Greengard work with banded matrices only. In Section 4, an 
important cancellation property of spectral integration will emerge from the explicit manner 
in which the linear systems are set up here. 

In this section, the interval of the boundary value problem is taken to be — 1 < y < 1 and 
the solution u is expanded as follows: 



The Chebyshev points yj = cos^jn/M) with j = . . . M are Af -|- 1 in number including the 
endpoints ±1, but the last coefficient in the Chebyshev series is always suppressed as indicated 
above. 

In each of the methods of this section and in the factored form of spectral integration 
described in the next section, the way the homogeneous solutions are computed may appear 
roundabout. As explained in Section 4, such a roundabout calculation is essential for ensuring 
accuracy. 



u{y) = Y ~^ "^'^^ ~^ ^ aM-iTu-i + O.Tm. 
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2.1 First and second order spectral integration 

Suppose a first order boundary value problem is given in the form 

{D - a)u = /. 



(2.2) 



The exact boundary condition is unimportant for much of the method. To begin with we 
assume the integral condition To{u) = or equivalently oq = 0. Suppose that the Chebyshev 
coefficients of / are given by fj = Tj{f). 



The indefinite integral of (2.2) gives 



u 



u + A 



/, 



(2.3) 



where A is an undetermined constant. By (|1.3|), the coefficient of T„ on the right hand side is 

Tn 



f 



fn-1 — fn+1 



2n 



for n= 1,2,...,M- 1. 



The n = M — 1 case assumes Jm = 0. Using (1.3) once again, the coefficient of r„ in the 



expansion of the left hand side of (2.3) is 



Tn u 



u + A 




2<n< M -2 

n = 1 

n = M -1 



for n = l,...,M — 1. The coefficients with n = 1 and n = M — 1 are obtained from the more 
general expression a„ — a(a„-i — an+i)/2n by setting oq = and om-i = 0, respectively. 
Equating coefficients for n = l,...,M — Iwe have M — 1 equations for the M — 1 unknowns 
ai, . . . ,aM~i- This tridiagonal linear system is solved to compute a particular solution 
satisfying {D — a)vP = f and the integral condition Tq{vP) = 0. 

A homogeneous solution v} satisfying {D — a)v} = and Tq{v}) = 1 is found as follows. 
We set n = 1/2 + u* so that u* satisfies To{u*) = and {D — a)u* = a/2. Thus u* is the 
particular solution of (2.2) satisfying 7o(ti*) = if / = —a/2 and it may be found using the 
method described for computing particular solutions. The same linear tridiagonal system is 
solved for computing u* and but with different right hand sides. 

The solution u is expressed as + Cu^ and the constant C is found using the boundary 
condition on u. 

Now we consider the second order problem 



(2.4) 



Integrating twice, we have 



u + b J u + c J J u + A + By = J J f. 

To find a particular solution, we assume the integral conditions 7o{u) = Ti{u) = or equiva- 

th( 

fn 



lently oq = ai = 0. By (1.4), the coefficient of T„ of the right hand side is 

fn-2 



Tn 



f 



fn+2 

4n(n-l) 2(n2-l) ' 4n(n + 1) 



+ 
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for n = 2, . . . , M — 1 (here /m = /m+i = is assumed). The coefficient of the left hand side is 
Tn(u + bfu + cf f u + A + By] = an-2 ( , , ^ — tt) + Qn-i ^ ° 



4n(n — 1) / \ln 



2(n2-l); "^'V2n; V4n(n + 1) 

for ?i = 2,...,M — 1. The vahdity of the equations for ?i = 2,3 reUes on the boundary 
conditions ai = 02 = 0. The equations for re = M — 2,M — 1 assume au = om+i = 0. 
The coefficients for n = 2,...,M — 1 on the left and right hand sides are equated to solve 
for the unknowns 02, • • • , olm-x- The particular solution uF obtained in this manner satisfies 
{jy^ + hD + c)vP = / and the integral conditions Tq{vP) = T\{yP) = 0. 

The first homogenous solution satisfies (D^ + hD + c)v^ = and the integral conditions 
Toiu^) = l,Ti{u^) = 0. To find it, we set = 1/2 -\- u*. Then u* satisfies the inhomogeneous 



equation (2.4) with / = — c/2 and the integral conditions 70(1**) = Ti{u*) = 0. The solution 
u* is computed using the same pentadiagonal system used for but with a different right 
hand side. 

The second homogenous solution satisfies the integral conditions To{v?) = 0,7i(u^) = 1. 
If we set v? = Ti + u*, u* satisfies the inhomogeneous equation with / = — (5 + cTi). It is 
found by solving the same pentadiagonal system. 

The solution u of (2.4) is expressed as + Cu^ + Dv? and the constants C and D are 



determined using the boundary conditions on u. 
2.2 Spectral integration of r-th order 

Define the operator L as Lu = u^"^^ +ar-iu^^~^^ + ■ ■ ■+aiu^^^ +aou. Consider the inhomogeneous 
equation Lu = f. A particular solution satisfying the integral conditions 

ro(n) = ... = %-i{u)=0 

or equivalently ao = ■ ■ ■ = a^-i = may be found as follows. The inhomogeneous equation is 
written in an integral form as 

/r^ r— 1 /.f- 

uH hao/ u + Y^ Ajy^ ^ / ■^^ 
j=o 

Using formulas for f'^Tn analogous to (1.3) and (1.4), we may express the coefficients of 



Tr, . . . , Tm-1 in terms of a.^, . . . , oaz-i. These coefficients are equated to the coefficients of 
/ and solved for a^, • . . ,aM-i to find the particular solution u^. The linear system has 
2r + 1 diagonals. 

To find the j-th homogeneous solution for j = 1, . . . , r, we first set = Tj^i + u*. The 
j-th homogenous solution satisfies the conditions Tk{u^) = 0, ifO<A;<r — 1 and A; 7^ j — 1, 
and Tj-i{u^) = 1. The function u* satisfies Lu* = —LTj and the first r coefficients in its 
Chebyshev series are zero. It can be found in the same manner as the particular solution. 

The solution of the linear boundary value problem is expressed as u^ + X]j=i CjuK The 
boundary conditions satisfied by u are used to determine the constants Cj. 
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Formulas for J"* T„ can be derived but get complicated. The 2r + 1 diagonal system can be 
difficult to set up correctly in programs. For the difficulties that arise for r = 4, see [9]. Spectral 
integration of r-th order will, however, prove quite useful in the discussion of cancellation errors 
in Section 4. 

2.3 Zebib-Greengard form of spectral integration 

The formulation of the Zebib-Greengard form of spectral integration given here makes it clear 
that the matrix systems solved have 2r + 1 diagonals if the problem is of order r. We assume 
L to be the operator defined in Section 2.2 above. 

The Zebib-Greengard form begins by assuming a Chebyshev series for u^''). We first find a 
particular solution of Lu = f subject to the integral conditions 

ro(n) = ro(n(i)) = • • • = ro(n(^-i)) = 0. 

When these conditions are used, the Chebyshev series of u^'^^ determines the Chebyshev series 
of u^^~^^ for s = r, . . . , 1 with no ambiguity. Normally, there is an undetermined constant of 
integration when the series of u^'^^ is integrated. But here the constant disappears because the 
mean mode of n^'^~^) is specified to be zero. Thus the Chebyshev series of Lu is determined 
unambiguously by the Chebyshev series of u^'^\ Coefficients in the Chebyshev series of Lu and 
/ are equated and solved for Tj (^u^'^^^ for j = 0, . . . , M. 

To find homogeneous solutions u, we expand u^'^^ in a Chebyshev series and take the 
integral conditions to be such that exactly one of To{u), . . . , To{u^^~^^) is one and the others 
are all zero. It is harder to find homogeneous solutions here than in the r-th order spectral 
integration method described in Section 2.2. One has to find polynomials pk of degree k for 
/c = 0, 1, . . . , r - 1 such that 7i)(pf ^) = 1 but Tkiplf^) = for d = 0, . . . /c - 1. 

This form of spectral integration has been derived in greater generality, allowing for linear 
operators with polynomial coefficients, in (TJ Section 3]. 

3 Spectral integration in factored form 

A linear operator L with constant and real coefficients can be factorized as 

L = {D-ai)...{D- am){D^ + hD + ci) . . . {D^ + bnD + Cn) 

where the coefficients are all real. We assume m + n > 2 and derive a method for solving 
Lu = f subject to boundary conditions that exploits this factorization of L. This method 
relies on spectral integration of orders one and two described in Section 2.1. The presentation 
of the method may appear more complicated but its implementation is much simpler than the 
methods of Sections 2.2 and 2.3. In addition, the factored form of spectral integration lends 
itself to piecewise Chebyshev grids as shown in Section 7. 

A particular solution is found by solving the following equations subject to integral condi- 



3 Spectral integration in factored form 



9 



tions on their solutions: 



{D - a2)u^ 



m+2n-2 



''m+2n-l 



7B(n^+2n-l) = 



{D - am)ul^ = u\^^^ %{u\^) = 



{p^ + hlD + Cl)uf„_2 = 7B('«2n-2) = ■n(w2n-2) = 



P 
U2 



roK)=riK) = o. 



This hst of equations is solved from first to last. Each equation is solved using one of the two 
methods described in Section 2.1. The subscripts on u, as in , indicate the number of 
"derivatives" in the function relative to Uq which satisfies Lu^ = f and is therefore a particular 
solution. 

If m > 1, the homogeneous solution with /i = 1 is found as follows. To begin with we 
solve the homogeneous problem 



{D - ai)u^+2m-i = To{u'^ 



+2m-l) 



1 



as described in Section 2.1. Thereafter, the inhomogeneous problems 



m-j 



u: 



n+2m—j+l 



To{u: 



'n+2m—j J 



are solved in the order j = 2, . . . ,m followed by the solution of 

(£)2 ^ ^ Ck)4n~2k = 4n-2k+2 To(4n-2fc) = 



n-2k) 



in the order k = 1, . . . ,n. The last solution to be found is = Uq and it satisfies Lu^ 



(3.1) 



(3.2) 



0. 



The inhomogeneous equations (3.1) and (3.2) are solved as described in Section 2.1. 

More generally, the homogeneous solution with 1 < h < m is solved beginning with the 
homogeneous problem 



2 \ -h 

ah)Um+2n-h 







To{u: 



'm+2n-h) 



1 



followed by the solution of (3.1 ) with j = /i + 1, . . . , m and (|3.2|) with k = 1, . . . ,n. As before, 
Uq is the last solution to be found and 



If h = m + 2i — 1 with 1 < i < n, the homogeneous problem solved at the beginning is 



{D^ + biD + Ci)u 



h 

2n-2i 







ro(4 



n-2i) 



1, ri(4n-2^) = 0. 



This is followed by the solution of (3.2 ) with k = i + . . . ,n. As before, Uq is the last solution 



to be found and = Uq. On the other hand, if h = m + 2i with 1 < i < n, the homogeneous 
problem solved at the beginning is 



{D^ + hD + Ci)u^2n-2r 

by the sc 
to be found and v!^ = Uq 



To{4n-2i) = 0, rMn-2^) = 1- 



This is followed by the solution of (3.2 ) with k = i + . . . ,n. As before, Uq is the last solution 
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By using the methods of Section 2.1 repeatedly, we end up with a particular solution 
and homogenous solutions u^, . . . , {t™"*"^". The solution of the boundary value problem Lu = f 
is expressed as 

m+2n 

u = uP + CjU^. 
i=i 

The constants Cj are found to fit the boundary conditions on u. 

There are two ways to find Cj. In the first method, the particular solution and the 
homogenous solutions are obtained in physical space as numerical values at the M + 1 points 
on the Chebyshev grid. Boundary conditions such as m(1) = A or u"(— 1) = B are expressed 
using a linear combinations of function values at the grid point. A boundary condition such as 
ii(l) = A simply specifies the function value at a single grid point. A boundary condition such 
as u"{—l) is interpreted as specifying that a certain linear combination of function values, the 
linear combination being determined by a single row of a spectral differentiation matrix [TI] , 
must have a specified value. This is the easier method for implementation and the one we 
have implemented. If the number M is not too large, this method will be adequate. If M is 
very large, then errors will creep in through the boundaries. In Section 7, we demonstrate that 
one can make do with small M even for computing solutions that develop very thin boundary 
layers or internal layers. 

The second technique uses the intermediate objects created when the particular solution 
and the homogeneous solutions are found. We illustrate the technique using an example. 
Suppose the boundary value problem is 

(D - ai){D - a2){D - a3){D - 04)^ = / 

subject to u{±l) = u'{±l) =0. If u = + J2'j=iCjU^ , the conditions on n(ibl) give two 
equations for the Cj after evaluation at ±1. We may rewrite the other boundary conditions as 
(D — 04)11 = at ±1. If we now note that 

3 

{D — 04)11 = Ui + ^ CjUi 
i=i 

we get two more equations for the Cj by evaluating at ±1. In the sum above, the j = ^ term 
does not appear. That is because the homogeneous solution satisfies {D — 04)^^ = 0. 

4 Errors in intermediate quantities and how they cancel 

Greengard |1] has emphasized that when spectral integration is used to solve a linear boundary 
value problem, the number of grid points should be adequate to resolve the solution but, most 
importantly, need not resolve the Green's function. For example, the solution of the linear 
boundary value problem {D^ — a?')u = — (vr^ + a?) sinvry with boundary conditions u(±l) = 
is n = sin Try. The solution can be represented with machine precision on a Chebyshev grid 
that uses slightly more than 20 points. If a = 10® it will take a Chebyshev grid with around 
10^ points to resolve the Green's function. Spectral integration can solve this boundary value 
problem using a Chebyshev grid with 20 or 30 points even if a = 10®. In this section, we explain 
the rather roundabout manner in which spectral integration comes to acquire this property. 
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If X is the solution of the matrix system Ax = 5i + 62 and xi, X2 are the solutions of 
Axi = bi, Ax2 = 62, then x = xi + X2 A is nonsingular. In machine arithmetic and in the 
presence of rounding errors, this linear superposition property will be true only approximately. 
This section deals with discretization errors and not rounding errors. Therefore we will assume 
this linear superposition property. 

Suppose that Lu = / is the given equation. With given boundary conditions on u, this 
equation is assumed to have a solution that is well-resolved using M + 1 Chebyshev points. 
In all forms of spectral integration, a particular solution satisfying Lu^ = / is found using 
some other global conditions on u^. Typically it will take many more points than M to resolve 
u"^. Thus the computed will be inaccurate. However, the approximation to u obtained by 
combining with the homogeneous solutions will retain its accuracy for reasons we will now 
explain. Before turning to the explanation, we mention that this situation, where a solution can 
be resolved more easily than the Green's function or a particular solution that satisfies integral 
conditions, occurs frequently in practice [Bl [3 [TJ] and it is important for an implementation 
to have this property. 

The explanation takes its simplest form for r-th order spectral integration described in 
Section 2.2 and it is with that method that we begin. Suppose L is a linear differential 
operator with constant coefficients and order r as in Section 2.2. We begin by denoting the 
computed solution of 

Lu = LTj 

with integral conditions To{u) = ••• = Tr-i{u) = by Uj for j = 0, ...,r — 1. Thus the 
Chebyshev series of Uj is obtained by solving a banded system with 2r + 1 diagonals and a 
right hand side that corresponds to the Chebyshev series of LTj integrated r times. The Uj 
will be typically quite inaccurate. We will show that the Uj occur in the particular solution 
and the homogeneous solutions in such a way that they cancel when an approximation to the 
solution of Lu = f with the given boundary conditions is computed. 

Let ue he the solution of Lue = f which satisfies the given boundary conditions and is 
accurate to machine precision with a Chebyshev series of M terms. We rewrite ue as 

= ^ + aiTi H h ar-lTr-l + UR 

where Tj{uji) = for j < r. We may rewrite / as 

/ = "^-^^0 H \- ar^iLTr-i + /r 

where /r = Lur. 

The particular solution of Lu = /r, which is computed by r-th order spectral integration 
is = ur. This is because ur satisfies the integral boundary conditions, the first r of 
its Chebyshev coefficients being zero, as well as Lu = Jr and can be represented to machine 
precision using a Chebyshev series of M terms. By linear superposition, the particular solution 
of Lu = f satisfying integral boundary conditions that is computed is given by 

On 

uP = —Uo + aiUi + ■■■ + ar-iUr-i + ur. (4.1) 

Homogeneous solutions of Lu = are computed such that Tj{u) = 1 but with the other 
r — 1 Chebyshev coefficients among the first r coefficients being zero, for j in the list 0, . . . , r — 
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1. This homogeneous solution is represented as u = Tj + u* and u* is computed as the 
particular solution of Lu* = —LTj, whose first r coefficients are zero. Therefore the computed 
homogeneous solutions are 

= l/2-Uo/2, u^ = Ti-Ui, u'' = Tr-i-Ur-i. (4.2) 



By observing (4.1) and ( |4.2[ ), we recognize that 

ri '^O —1 —2 — r 

Ue = + —U + aiU + • • • + Or-lU . 



In this linear combination of the particular solution with the homogeneous solutions, the 
coefficients are such that the inaccurate Uj cancel exactly and the solution ue satisfies the given 
boundary conditions. If the equations that are solved to determine the linear combination of 
homogeneous solutions with the particular solution are reasonably well-conditioned, which we 
may expect because these are typically very small linear systems, the computed solution will 
produce ue very accurately. 

The explanations for the factored form of spectral integration and the Zebib-Greengard 
version are more complicated. We will give the explanation for the problem {D—a){D—b)u = f. 
The given boundary conditions are assumed to be ^(±1) = 0. We assume as before that ue 
is the approximate solution whose Chebyshev series has M terms and which is accurate to 
machine precision. 

Suppose U2 is the solution oi {D — b)u = satisfying 7o(f^2) = computed as explained in 
Section 2.1 using a Chebyshev series with M terms. Similarly, let U[ be the computed solution 
of (D — a)u = 0, and let U2 be the particular solution of {D — h)u = U[ satisfying To{U2) = 
and computed using a Chebyshev series with M terms only. The factored form of spectral 
integration uses Ui and U2 as a basis of homogeneous solutions. For reasons given above, Ui 
and U2 are typically very inaccurate. 

As before, we will split ue but the split is more complicated this time. We write 

UE = Y + i'^l ~ 7)7"! + UR, 

where 7 is chosen such that 

UR = jTi + a2T2 H h um-iTm-i 

satisfies {D — b)uR = 0. By applying {D — a){D — b) to ue, f can be split as 

/ = ^ - (a + 6)(ai - 7) + ab{ai - ^)Ti + {D - a){D - b)uR. 

The computed particular solution that corresponds to {D—a)(D—b)u = abao/2 — {a+b){ai —7) 
is {abao/2 — (a + b){ai — 7)) C/i. The particular solution of 

{D-a){D-b)u = ab{ai-j)Ti (4.3) 

is obtained by solving 

{D - a)v = ab{ai - -f)Ti = {D - a) {-b{ai - j)Ti) + b{ai - 7) 
{D - b)u = V. 
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Because of the way the right hand side of the {D — a)v equation is rewritten, the particular 



solution of (4.3) may be taken to be computed as the particular solution of 



{D - h)u = -h{ai - 7)ri + b{ai - ^)U[ = {D - b){ai - j)Ti - (ai - 7) + 5(ai - -f)U[. 



From the form of the right hand side, we infer that the particular solution of (4.3 ) is computed 
to be 

(ai - 7)ri - (ai - 7)C/2 + K«i - 
Because of the way / was split, 

uP = {abao/2 - (a + b){ai - 7)) Ui + {ai - j)Ti - {ai - ^)U2 + b{ai - j)Ui + ur 

= a {bao/2 - ai + 7) C/i - {ai - j)U2 + (ai - j)Ti + ur (4.4) 

is the particular solution of {D — a){D — b)u = f computed by the factored form of spectral 
integration. 

The homogeneous solutions computed by the factored form of spectral integration are 

-1 aUi , U2 

^=2 + ^- 



By observing (4.4) and (4.5), we find that 

ue = vP — 2 {bao/2 — ai + 7) 'U"'^ + aou^. 

We may argue as before that even though n^, u^, v? are inaccurate, the factored form of 
spectral integration solves {D — a){D — b)u = f with boundary conditions u{H) = accurately. 

5 Condition numbers and numerical accuracy 

Spectral integration of r-th degree solves a banded matrix with 2r + 1 diagonals. The same 
is true for the Greengard-Zebib version. In its factored form, the matrices solved by spectral 
integration have either 3 or 5 diagonals. In all these cases, the condition number of the matrix 
remains bounded asM— )■ 00 [Hll], M + 1 being the number points in the Chebyshev grid. The 



matrices discretize a compact integral operator such as (1.8). Therefore the condition number 
does not diverge as the number of grid points is increased. 

A question of greater practical import is the limit the condition number converges to as 
M —7- 00. Numerical tests show that for an operator such as (D'^—a^) the limit is approximately 

for large a (see the last columns of Tables 2 and 3 of [1 ). For an operator such as (D^ — 
a^)(Z)^ — 6^) the limit appears to be a^6^. If a = 10^ and 6 = 2 x 10^, for example, the condition 
numbers will be greater than 10^^. Yet the boundary value problem [D"^ — o?){D'^ — b'^)u = o^b"^, 
li(ibl) = li'(ibl) = 0, is solved with a loss of around 6 digits. The boundary layer has a thickness 
of approximately 10~^ and therefore a loss of around 6 digits is inevitable. The surprising fact 
is that all versions of spectral integration are much more accurate than the condition number 
suggests. 



Figure 5.1 shows the singular values and vectors of a spectral integration matrix. It may 



be observed that the first few singular vectors (the right singular vectors are shown but plots 
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Fig. 5.1: Singular values and vectors of a spectral integration matrix corresponding to + 
aD + h with a = 10^, h = —10^, and M = 128. The plots of the singular vectors show 
absolute value of the entry against its position in the vector. 
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a 


M 


error 




a 


M 


error 


le+06 


1024 


0.497459 




le+06 


8 


6.15934e-09 


le+06 


4096 


7.22283e-05 




le+06 


16 


4.5C-16 


le+06 


8192 


3.81267e-ll 




le+06 


32 


4.11129C-16 


le+06 


16384 


4.9402e-ll 




le+06 


1024 


5.98272e-14 


le+06 


65536 


1.55942e-10 




le+06 


65536 


2.3926e-12 


(a) {D + a)u = a with m(-1) = 0. 


(b) {D + a 
with u{ 


)u = -Kcos{'Ky) + osin(7rj/' 
-1) = 0. 


a 


M 


error 




c 


M 


error 


le+06 


1024 


0.391632 




10000 


8 


4.21405e-06 


le+06 


4096 


9.28769e-05 




10000 


16 


2.27521e-12 


le+06 


8192 


4.81313e-ll 




10000 


32 


7.44196e-15 


le+06 


16384 


4.63216e-ll 




10000 


1024 


2.60902e-14 


le+06 


32768 


4.25722e-ll 




10000 


16384 


2.74558e-13 



(c) {D^ - aD)u = with u(-l) = 1 and (d) {D^ + c)u = (-tt^ + c) sin ttj/ with 

u{l) = 2. m(±1) = 0. 



a 


b 


M 


error 1 


error2 


le+06 


2e+06 


1024 


0.863351 


0.863351 


le+06 


2e+06 


8192 


2.14342C-07 


2.14697e-07 


le+06 


2e+06 


16384 


1.11927e-09 


8.68444e-10 


le+06 


2e+06 


131072 


2.62727e-08 


3.47769e-08 



(e) - a'){D' - b'')u = with m(±1) = u'{±\) = 0. 



Tab. 1: Errors in spectral integration. The version of spectral integration used is the factored 
form of Section 3. In (e), the two errors correspond to the factorizations {D — a){D + 
a){D - h){D + h) and (D^ - a'^){D'^ - 6^), respectively. 

of the left singular vectors look similar) are localized with most of their energy in the first few 
entries. In all versions of spectral integration, the right hand side vector has entries whose 
magnitudes decrease rapidly after the first few entries. Therefore the localization of singular 
vectors may be the reason solutions which develop boundary layers are much more accurate 
than the condition numbers of linear systems that occur in spectral integration suggest. 

6 Numerical examples 

Spectral integration can solve linear boundary value problems using no more than the number 
of grid points sufficient for resolving the solution. For reasons explained in Section 4, it is not 
necessary for the grid to resolve the Green's function of the linear boundary value problem. 
However, if the grid over-resolves the solution, spectral integration remains accurate. In higher 
order problems whose solutions develop thin boundary layers, spectral integration is far more 
accurate than indicated by the condition numbers of matrices that are solved. This could be 
because of the phenomenon of localization of singular vectors discussed in the previous section. 
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Table [T] applies the version of spectral integration derived in Section 3 to some examples and 
demonstrates these properties. 

Example (a) of Table [T] develops a thin boundary layer at y = —1. The boundary layer 
thickness is approximately 10~^. Within this thin boundary layer, the solution transitions from 
a value approximately equal to 1 in the interior of the domain to the value of at y = —1. The 
table shows that the solution is well resolved with M = 8192, where M + 1 is the number of 
points in the Chebyshev grid. If the grid were uniformly spaced, it would take about a million 
points before the grid saw the boundary layer. Thanks to quadratic clustering of Chebyshev 
points near the end points ±1, far fewer points are used. If the boundary layer is as thin as 
10~^, any numerical method will loose about 6 digits of accuracy. The boundary layer will 
look like 1 — e""*-^"*"^^ near y ~ —1. The mere evaluation of 1 — e~°^^^^^ near y = —1 will loose 
about 6 digits if a = 10®. Thus the errors in Table [l^ are the best possible in double precision 
arithmetic. All errors reported here are in the L°° or sup norm. 

Example (b) of Table [T] uses the same operator as example (a). However, the right hand 
side is chosen such that the solution is sinvry. We find that the solution is computed with 
accuracy comparable to machine precision with M = 16. Thus the grid is required to resolve 
the solution but need not resolve the Green's function. The reason spectral integration has 
this important property is discussed thoroughly in Section 4. Examples (a) and (b) as well 
as all the other examples demonstrate the robustness of spectral integration when the grid 
over-resolves the solution. In practice, this property of spectral integration is not as important 
as the others. 

Example (c) is taken from [5] . In the next section. Table [Tj^c) will be compared to spectral 
integration using a piecewise Chebyshev grid. We will find that the piecewise Chebyshev grid 
has dramatic advantages. 

Example (d) demonstrates that even when the Green's function has rapid oscillations, the 
grid has to resolve the solution but not the Green's function. 

Example (e) treats a fourth order boundary value problem. The two errors reported in 
Table[l]3 correspond to the factorizations {D - a){D + a){D - b){D + b) and {D'^ - a?){D'^ 
either of which can be used with the method of Section 3. There is no significant difference in 
errors between the two factorizations. 

Table [2] compares a Chebfun implementation of spectral integration with C-|— 1-. The Cheb- 
fun system supports functions as first class objects in MATLAB [13]. The Zebib-Greengard 
version of spectral integration has been implemented in MATLAB using chebfuns P]. The 
Chebfun implementation does not use Chebyshev series as in Sections 2 and 3 but works in 
physical space using integration matrices that are presumably based on the Clenshaw-Curtis 
rule. The Chebfun implementation was unable to solve the examples (c) and (e) of Table [T| It 
detects that the matrices are ill-conditioned and aborts the calculation. The first two rows of 
Table |2] show that the Chebfun errors are noticeably greater than the C-|— |- errors. It appears 
that the partial immunity to ill-conditioning of matrices enjoyed by the methods of Sections 2 
and 3 is not possessed by the Chebfun implementation. That is probably because the singular 
vectors are localized when the matrices solve for Chebyshev coefficients, as we illustrated in 
Section 5, but are not localized in physical space. 

The last row of Table |2]shows that C-|--|- is 15, 000 or 60, 000 times faster than the Chebfun 
implementation. Before we interpret those numbers, we point out that Chebfun is a general 
purpose system than can solve a wide variety of problems. It is not designed to be fast. 
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Problem 


Chebfun 


C++ 


(D'^ - aD)u = with 


Error of 2.68 x 10"*^ with 


Error of 2.47 x 10"^^ with 


u(-l) = 1, n(l) = 2, and 


M = 744. 


M = 744. 


a = 10"^. 






{D'^ + 5D + 10^)^ = 


Error of 1.71 x 10"^^ with 


Error of 1.50 x 10^^^ with 


-500cos(100y)e~^^ with 


M = 99. 


M = 100. 


n(0) = and 






= sin(100)e-^ 






{D'^ - a^){D^ - b'^)u = f 


1.5523 billion cycles for a 


108,600 or 25,900 cycles 


with / and the boundary 


single solve on a single 


for a single solve on a 


conditions chosen such 


core of 2.67 GHz Xeon 


single core of 2.67 GHz 


that u = sinvry. The 


5650. 


Xeon 5650. 


parameters a and b are 






0{1). 







Tab. 2: Comparison of DriscoU's Chebfun/MATLAB implementation of spectral integration 
with a C++ implementation of the method in Section 3. 



Nevertheless, discussion of why the speed-up in C++ is so large is pertinent because credible 
information about the speed of MATLAB with respect to C++ is so hard to find and the 
folklore beliefs on this topic are spurious. 

Programs written in C++ cannot be benchmarked perfectly. The execution time of the 
same instruction sequence depends upon contextual and architectural factors such as cache 
misses, the load on internal buffers used by the memory controller, the load on reservation 
stations used for register renaming, the efficacy of DRAM optimizations, and other factors. 
The manner in which reliable timing numbers can be obtained is discussed briefly in 
Reasonable precautions have been taken to ensure the validity of numbers reported in Table 

El 

In [TT], we compared carefully optimized C++ programs with MATLAB and found C++ 
to be a few thousand times faster. The comparison given in Table [2] is perhaps more realistic. 
The C++ program was written in the ordinary way, with no effort to optimize by unrolling 
loops and inspecting the assembly code as in [11_ , although we used good libraries following our 
usual practice. We did not decouple even and odd modes, an optimization that the problem 
in row three of Table |2] permits. The MATLAB program, which was kindly supplied to us by 
DriscoU [2;, is expertly coded using an idiom that is native to that system. 

We now return to the C++ speed-up factors relative to Chebfun/MATLAB, which are 
implied to be 15,000 and 60,000 by Table |2] . The speed-up of 60,000 comes about if the 
differential operator is fixed but the right hand side is varied. Our C++ program is written to 
be optimal in that situation. However, in Chebfun/MATLAB each solution of the boundary 
value problem builds the matrices from scratch. The speed-up factor of 15,000 was obtained 
when we did the same. Specifically, the parameters in the differential operator were varied 
randomly from one use to the next which meant that matrices needed to be set up anew and 
factored for each usage. 

The Chebfun/MATLAB algorithm is different from the algorithm implemented in C++. 
In particular, Chebfun tries to solve a harder problem by attempting to figure out the right 



7 Piecewise Chebyshev grids 



18 



number of grid points to resolve the solution on the fly. That information is supplied to the 
C++ program explicitly. In [lO], it is reported that Chebfun computes a representation of 
(.01 + (x - 0.3)2)"^ in 0.1 seconds on a desktop. We take that to be about 0.25 billion cycles. 
The number of cycles expended by Chebfun to determine the collocation points automatically 
for the test problem in row 3 of Table |2] is certainly smaller than 0.25 billion cycles. The 
right hand side in the test problem is of the form A + B sin ny + C cos iry and the solution is 
sin ny, both of which are easy to resolve. Chebfun's ability to collocate automatically is not a 
significant influence on the observed speed-up. 

The resolution used by Chebfun in the third row of Table [2] was M = 22. We used 
M = 32 in C++. Chebfun used dense matrices while our C++ program uses the discrete cosine 
transform (DCT) and banded matrices. The real-to-real DCT has a substantial constant factor 
for M = 32, and our opinion is that this difference accounts for a speed-up factor of at most 10. 
We also mention that using banded matrices and the DCT is not native to MATLAB's idiom. 
While such optimizations are natural in C++, their use would make MATLAB programming 
cumbersome and error prone. 

Thus we are led to the conclusion that the performance penalty incurred by the use of 
MATLAB is a factor of a few thousand or more. A similar performance penalty was observed 
in [llj. This performance penalty is likely to increase in the next few years and we expect it 
to be much greater for parallel programs. Progress in computer architecture has far out-paced 
progress in software. A C++ program written without awareness of computer architecture 
had a much better chance of approaching optimality twenty years ago than it does today. The 
import of this trend can be observed within MATLAB itself. In m] , we reported that a MAT- 
LAB function for computing Lagrange weights was six times faster when written with explicit 
loops, in a style that is natural to C or C++ but not to MATLAB, instead of vectorization, 
which is the idiom native to MATLAB. Presumably MATLAB uses some kind of a byte code 
interpreter and explicit loops convert to better byte code than vectorized statements. 

7 Piecewise Chebyshev grids 

Section 7.1 generalizes spectral integration to piecewise Chebyshev grids. The interval of the 
boundary value problem, which we usually take to be [—1,1], is split into sub-intervals and 
each sub-interval is discretized using scaled and shifted Chebyshev points. We find that a 
problem that requires M = 8192 if the entire interval is discretized using Chebyshev points 
can be solved with the same level of accuracy using three intervals and Mi = M2 = = 32. 

If M is as small as 32, differentiation matrices are as accurate as spectral integration. 
They are much simpler to implement although they use more memory and more arithmetic 
operations. In Section 7.2, we show that boundary layers and internal layers can be computed 
accurately using differentiation matrices and piecewise Chebyshev grids. 

7.1 Spectral integration and piecewise Chebyshev grids 

To generalize spectral integration to piecewise Chebyshev grids, we consider the operator 
{D — a){D — b) over the interval — 1 < y < 1 and the boundary value problem {D — a){D — h) f = 
u. As in Sections 2 and 3, the boundary conditions enter only at the end and much of the 
method is independent of the specific form of the boundary conditions. The generalization to 
operators of the form {D — ai) . . . {D — a^) will be obvious. 
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Let [—1,1] = Xi U . . . U X^, where Xj are intervals with disjoint interiors and with the 
right end point of Xj equal to the left end point of Xj+i. Thus Xi, . . . ,X„ are disjoint intervals 
arranged in order. Let Wi denote the width of the interval Xj. 

We use a linear change of variables Xj — )■ [—1,1] and rewrite the given differential equation 

as 

f bwi\ wf 



after the change of variables. In (7.1 ) it is assumed that u and / have been shifted from Xj to 
[—1,1] although that is not indicated explicitly by the notation. 
We define Ui as 

Ui = uP + ax,u^ + (7.2) 

where is the particular solution and , v? are the homogeneous solutions of (7.1), computed 
as described in Section 3. 

For i = 1, . . . ,n, the coefficients aj. and comprise 2n unknown variables in total. We 
will solve for these unknowns using the two boundary conditions and continuity conditions 
between intervals. The boundary conditions give two equations such as 

ni(— 1) = left value 
n„(l) = right value. 

For i = l,...,n — 1, the continuity conditions are 

Ui{l) = Uj+i(-l) 

{{D-Wib/2)ui){l) ^ {{D - w^+ih/2)ui+i) (-1) 

The second continuity condition requires the derivatives to be continuous while accounting for 
the shifting and scaling of intervals of width Wi and Wi+i to [—1, 1]. The function [D — Wih/2)ui 
is available through the intermediate quantities generated by the method of Section 3. In 
particular, we have, 

{D - Wib/2)u'P = u\, {D -Wih/2)u^ = u\, {D - Wib/2)u'^ = 



in interval Xj. Once we solve for aj. and for i = 1, . . . , n, we may use (7.2) to form Uj. The 
solution u is obtained by shifting the Ui from [—1,1] back to Xj. 

Table |3] summarizes spectral integration of piecewise Chebyshev grids applied to solve 
(D^ — aD)u = with n(— 1) = 1, u{l) = 2, and a = 10^. The boundary layer of this problem 
was displayed in Figure |1.1[ It is evident from the table, that the intervals must be chosen 
carefully. We were unable to get an accurate solution with fewer than a thousand grid points 
and just a single interval properly contained in the boundary layer. The last row of Table [3] 
reports a solution with Mi = M2 = M3 = 32 and an error of 4.7 x 10~^^. In that computation, 
two intervals are contained inside the boundary layer. In Table [l]: we found that a single 
Chebyshev grid required M = 8192 to reach more than ten digits of accuracy. 

The example of Table |3] coincides with Example 3 of Table 7 of that important paper 
reports an error of 2.33 x 10~^^ using 15 intervals and M = 16 for each interval (the total 
number of grid points is 321). As discussed in the introduction, the method of [5, can handle 
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Ml 


M2 


M3 


node 2 


node3 


error 1 


error2 


16 


1024 


32 


0.5 


0.99999 


5.80845e-06 


1.701397e-02 


16 


4096 


32 


0.5 


0.99999 


4.07361e-ll 


2.825425e-10 


32 


128 


32 


0.999 


0.99999 


4.49718e-ll 


9.037897e-10 


32 


64 


32 


0.9999 


0.99999 


4.33247e-ll 


1.349214e-10 


32 


32 


32 


0.99995 


0.99999 


4.66069e-ll 


8.507683e-ll 



Tab. 3: Solution of — aD = with u(— 1) = 1, u(l) = 2, and a = 10^ using a grid with 
three intervals, which are discretized using Ml + 1, M2 + 1, and M3 + 1 Chebyshev 
points, respectively. Nodes 1 and 4 are located at ±1. Errors 1 and 2 correspond to 
spectral integration and the use of differentiation matrices, respectively. 




-8 -6 -4 -2 2 4 6 8 



(a) (b) 

Fig. 7.1: Solution of u" + ayu' = with u{—l) = —1, n(l) = 1, and a = 10^'^ using dif- 
ferentiation matrices and a piecewise Chebyshev grid, (a): Spy plot of matrix with 
Ml, M2, M3, M4, M5 = 4, 6, 8, 10, 12. (b): Transition region of the solution. 

non-constant coefficients but is significantly more expensive than spectral integration when the 
coefficients are constant. 

Shiskin meshes [8J and finite element discretizations can be used to resolve boundary layers. 
For example, the boundary value problem —e^u" + u = y with boundary conditions it(±l) = 
and e = 10~^ has been solved using a Shiskin mesh by Zhang jTT]. An error of 2.8 x 10"^ 
was realized using 1024 degrees of freedom and linear elements. An error of 3.14 x 10~^ was 
realized using 512 degrees of freedom and quadratic elements. The thickness of the boundary 
layers is approximately e or 10~^ for this example. Shiskin meshes approximately locate the 
transition points from the internal region to the boundary layer and make the mesh finer and 
finer as the boundary is approached. An apriori bound on L°° error at some but not all of the 
grid points has been obtained [SI EE] • 
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7.2 DifFerentiation matrices and piecewise Chebyshev grids 

The final item is the use of piecewise Chebyshev grids and spectral differentiation matrices. 
The ease with which a variety of problems can be solved using spectral differentiation matrices 
is illustrated in An efficient and stable method for computing spectral differentiation 

matrices is given in . Spectral differentiation matrices become impractical and highly inac- 
curate in the presence of boundary layers or internal layers. However, if piecewise Chebyshev 
grids are used the number of grid points on each sub-interval can be kept small and the use of 
spectral differentiation matrices becomes possible. We illustrate using two examples. 

The first example is the same as in Table [sj namely, {D^ - aD) u = with u{-l) = 1, 
u{l) = 2, and a = 10^. We again use three intervals Ti = [—1, 1^2] , I2 = [??2, "'^s] , and X3 = [7^3, 1] . 
Each interval is discretized by scaling and shifting the Chebyshev points. We assume the 
number of Chebyshev points in the intervals to be Mi -|- 1, M2 + 1, and M3 -|- 1, respectively. 
Since the boundary conditions specify u at the end points, and u must be continuous at r/2, r/3, 
the total number of unknowns is Mi -|- M2 -|- M3 — 1. We denote the vector of values of the 
function u at grid points by U . Its entries are the unknown values of u at interior grid points 
and the two known values at ?/ = ±1. 

Denote the differentiation matrix over the interval Xj by and the matrix corresponding 
to (Z?^ — aD) by £j. Both Pj and Ci are matrices of dimension (Mj -|- 1) x (Mj -|- 1). 

Over Xi, the given equation can be written as CiUxi = 0, where Uxi is a subvector of U 
corresponding to the Mi -|- 1 grid points in Xi. Only Mi of the entries of C/j^ are unknown 
because the value of ti at —1, which is the left end point of Xi, is known to be 1. We enforce 
this equation at Mi — 1 interior points of Xi only. In particular, of the M-|- 1 equations CiUx^ 
the first and last are dropped. This gives us Mi — 1 equations. 

Similarly, the equations C2UX2 = and C^Ux^ = are enforced at the M2 — 1 and M3 — 1 
interior points of Xi and X2 only. In the case of C^Ux^ = we must take note of the fact that 
the value of n at 1, which is the right end point of X3 is already known. 

In this way, we get Mi -|- M2 -|- M3 — 3 equations for Mi -|- M2 -|- M3 — 1 unknowns. One more 
equation is obtained by setting the entries of ViUxi and V2UX2 that correspond to the derivative 
at r/2 equal to each other. Yet another equation is obtained by matching the derivatives at r/3. 

Table |3] shows that solutions obtained using differentiation matrices and a piecewise grid 
are as accurate as the solutions from spectral integration. 

The second example we consider is eu" + yu' = with boundary conditions ^(±1) = ±1 
and e = 10~^^. The exact solution of this boundary value problem is given by 



u{y) = -1 + 



2 fy/^' e-*' dt 



The solution has an internal layer at y = of width approximately e^^/^ or 10^®. In Figure 



7.1 , we show the spy plot of a matrix corresponding to division of [—1, 1] into five sub-intervals 
as well as the transition region of the solution. 

This second example occurs near the end of [T] , where it is reported that mapped Chebyshev 
points with M = 1024 compute the solution with an overshoot of 3 x 10"'^. From Table |4] we 
see that the overshoot is reduced to the order of machine precision using only 161 grid points. 
The overshoot is seen to be highly sensitive to the location of the nodes. The solution plotted 



8 Acknowledgments 



22 



m 


node4 


overshoot 


32 


5^ 


3.7e-15 


32 


3Vi 


1.2e-08 


32 




8.6e-09 


24 


5^^ 


1.8e-08 



Tab. 4: Solution oi e u" + y u' = with u(— 1) = —1, u(l) = 1, and e = 10 is shown in Table 



7.1 



This table gives the overshoot beyond [—1,1] of the solution computed using 6 
nodes and 5 intervals. Nodes 1, 2, 3, 5, and 6 are fixed at —1, —8e^^'^, — 3e^/^, 8e^^'^, 
and 1, respectively. The number of Chebyshev points in each interval is m+ 1. 



in Figure |7.1| 3 corresponds to the top row of the table. The solution appears to have around 
10 digits of accuracy. The sensitivity to the placement of nodes can probably be reduced using 
Green's functions, which is the approach advocated in [5] and [T5] . 

To summarize, this paper has introduced new formulations of spectral integration that are 
particularly convenient for higher order problems. We have illustrated the partial immunity 
of spectral integration to ill-conditioning of linear systems that are solved. We have clarified 
the manner in which spectral integration produces an accurate final answer even when the 
intermediate quantities are inaccurate. Finally, we have illustrated the efficacy of piecewise 
Chebyshev grids for resolving boundary layers and internal layers. 
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